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In this paper, we discuss the algorithms used in the LO evolution program 
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^ . I. INTRODUCTION 

Due to the recent availability of exclusive hard diffraction data at HERA, there has been 
a great interest in the study of generalized parton distributions also known as nondiagonal, 
off-forward or non- forward parton distributions occurring in these reactions (see Ref. [1-11]). 
These parton distributions are different from the usual, diagonal distributions found in e.g. 
inclusive DIS since one has a finite momentum transfer to the proton due to the exclusive 
nature of the reactions. In this paper we give an exposition of the algorithms used to nu- 
merically solve the generalized GLAP-evolution equations. The main part of the evolution 
program was taken over from the CTEQ package for the diagonal parton distributions from 
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inclusive reactions. At this point in time the evolution kernels for generalized parton distri- 
butions are known only to leading order in a s and thus our analysis will be a leading order 
one. 

The paper is organized in the following way. In Sec. II we will quickly review the formal 
expressions for the parton distributions and the evolution equations together with the explicit 
expressions for the kernels and a first comment on the arising numerical problems. In Sec. 

III we will explain the difference of our algorithms to the ones used in the original CTEQ 
package and then give a detailed account of how we implemented our algorithms. In Sec. 

IV we demonstrate the stability of our code and show that we reproduce the case of the 
usual or diagonal parton distributions within 1% for a vanishing asymmetry factor. Sec. V 
contains concluding remarks. 

II. REVIEW OF NONDI AGONAL PARTON DISTRIBUTIONS, EVOLUTION 

EQUATIONS AND KERNELS 

A. Nondiagonal Parton Distributions 

Generalized or, from now on nondiagonal parton distributions, occur for example in 
exclusive, hard diffractive J/ip or p meson production and alternatively in deeply virtual 
Compton scattering (DVCS), where a real photon is produced. As mentioned in Sec. I since 
one imposes the condition of exclusiveness on top of the diffraction condition, one has a 
kinematic situation in which there is a non-zero momentum transfer onto the target proton 
as evidenced by the lowest order "handbag" diagram of DVCS in Fig. 1. The picture serves 
to only introduce the kinematic notations used throughout the text and nothing more. For 
more on DVCS see for example Ref. [6,7,12-20]. 

The nondiagonal quark and gluon distributions have the following formal definition as 
matrix elements of bilocal, path-ordered and renormalized quark and gluon operators sand- 
wiched between different momentum states of the proton as in the factorization theorems 
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FIG. 1. The lowest order handbag contribution to DVCS with Q 2 = —q 2 and q' 2 = 0. 
for exclusive vector meson production [4] and DVCS [7,19,20]: 

fq/p = r^ e ' lX2P+y '(p\ T ^y^^h +p ^\p')^ 

fg/p = ~ ^-^-^e- lX2P+y ~(p\TGt(0,y-,0 ± )PG»+(0)\p'}. (1) 

J — oo Z71 X1X2P 

with X2 = x\ — A where the asymmetry or nondiagonality parameter A is usually xbj in, 
for example, DVCS or exclusive vector meson production however not in diffractive di-muon 
production. 



B. The GLAP-Evolution Equations for Nondiagonal Parton Distributions 

The GLAP-evolution equations follow from the usual renormalization group invariance 
of the factorization formula and lead to the following evolution equations for the singlet (S) 
and non-singlet (NS) case [5-7]: 

dq N s(xi, A,Q 2 ) f 1 dy x 2 
d^Q 2 = L ^ qNS ^ Q ^ 

d9si 2^ Q2) = L [ p ^(m, A, Ql) + P 9 Mvu A, Ql)} , 
dqs{ 2^ Q2) = £ [ p ^s(m, a, Ql) + P q99s {Vu a, Ql)] . (2) 

Note that qs,NS = xifq/p, 9s = and the kernels to leading order [21] are given by 

[5-7]: 



3 



Pqq,S,Ns(Xl,&) 
P. 



C f 



IT 



a. 



qg,S{Xi, A) = — N F 
7T 



( Xl ,A) 



7T 



x\ + x\ — A(xi + x 2 ) 

(1-A)(1-X!) 

\x\ + xi(l - xi) 2 - x 2 A] 
(1 - A) 2 ' 
[i + (i-m) 2 -A] 
1- A 

;i-^i) 2 + (i-^)(^i-A) 



-(J(l-Xi) 



r 1 d^i Z" 1 dz 2 3 
/o z\ Jo z 2 2 



a) = n c [2 ^ ~ iy ;; 2 A ;; i/v ^ — ^ - 1 - ^ + , + 



xi - A 



(1-A) 2 
2iV c io 



dzi dz 2 

JO Zi io Zo 



1-xi (l-xO(l-A) 

(3) 



With our definitions, we obtain for the diagonal limit, i.e. , A = 0, qs,NS — > xQ(x, Q 2 ) and 
<7s — > xG(x, Q 2 ) where Q and G are the usual parton densities. 

A word concerning the above employed regularization prescription which is the usual + 
- prescription in the first integral below and a generalized + - prescription for the second 
integral, is in order, since these prescriptions have direct implications on the numerical 
treatment of the integrals involved. In convoluting the above kernels, after appropriate 
scaling of X\ and A with y 1 , with a nondiagonal parton density, one has to replace Z\ and 
z 2 in the regularization integrals with z\ — > (yi — Xi)/y 1 and z 2 — > (yi — x 1 )/(y 1 — A). This 
leads to the following regularization prescription as employed in our modified version of the 
CTEQ package and in agreement with Ref. [7]: 

1 dyx f(yi) f l dy l yif{yi) - x l f{x l ) 



L 



xi j/i 1 - x 1 /y 1+ 
d (x 1 -A)f(y 1 ) 
*i 1 {yi -x 1 )(y 1 - A) 



+ f(x) ln(l-X!) 
xi V\ Vi - xi 

dyi Vif(yi) - x 1 f(x 1 ) r l dy x y x j{v\) - A/(xi) 



(4) 



J x-i 
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Vi ~ xi 
1 — X\ 
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xi yi 



(5) 



Eq. 5 and a closer inspection of Eq. 3 reveals that if one were to integrate each term 
by itself one would encounter infinities in all the expressions at both the lower bound of 
integration if A = y 1 and in taking the limit A = x\. Although Eq. 3 is completely 
analytical, it will cause numerical problems since the cancellations of the infinite terms can 
only be done in the analytical expressions. This is in contrast to the diagonal case where 
such problems are absent. The integration over Q 2 is identical to the diagonal case and 



hence has already been dealt with in the original CTEQ-code. 

III. DIFFERENCES BETWEEN THE CTEQ AND OUR ALGORITHMS 

Let us point out in the beginning that our code is to 99% the original CTEQ-code (for 
a detailed account of this code see Ref. [22]). We only modified the subroutines NSRHSM, 
NSRHSP and SNRHS within the subroutine EVOLVE and added the subroutines NEWAR- 
RAY and NINTEGR. These routines are only dealing with the convolution integrals but 
not with, for example, the Q 2 -integration or any other part of the CTEQ-code which re- 
mains unchanged. This is due to the fact that the main difference between the diagonal and 
nondiagonal evolution stems from the different kernels which only influence the convolution 
integration and nothing else. 

In order to make the simple changes in the existing routines more obvious we will first 
deal with the new subroutines. 

A. NEWARRAY and NINTEGR 

Due to the increased complexity of the convolution integrals as compared to the diagonal 
case as pointed out in Sec. II B, we were forced to slightly change the very elegant and fast 
integration routines employed in the original CTEQ-code. The basic idea, very close to the 
one in the CTEQ-code, is the following: Within the CTEQ package, the parton distributions 
are given on a dynamical x- and Q-grid of variable size where the convolution of the kernels 
with the initial distribution is performed on the x-grid. Due to the possibility of singular 
behavior of the integrands, we perform the convolution integrals by first splitting up the 
region of integration according to the number of grid points in x, analytically integrating 
between two grid points X; L and x i+ i where i runs from 1 to the specified number of points 
in x and then adding up the contributions from the small intervals as exemplified in the 
following equation: 
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f ^lf( Xl / yi ,A/y 1 ,y 1 ) = Sf =0 f" 1 —f(x 1 /y 1 ,A/y 1 ,y 1 ), (6) 
J xi yi Jx, yi 

where f(x\/yi,A/yi,yi) is the product of the initial distribution for each evolution step and 
an evolution kernel with x = x±, x^ = 1. We can do the integration analytically between 
two neighboring grid points by approximating the distribution function f(yi) through a 
second order polynomial ay\ + by\ + c, using the fact that we know the function on the grid 
points Xi-i, Xi and x i+ i and can thus compute the coefficients a,b,c of the polynomial in the 
following way, given the function is well behaved and the neighboring grid points are close 
together [23]: 

f(x 1+1 ) = ax 2 i+1 + bx i+1 + c 

f(xi) = ax\ + bxi + c 
f(xi-i) = + bxi-i + c (7) 

which yields a 3 x 3 matrix relating the coefficients of the polynomial to the values of 
the distribution functions at Xi_i,Xi and x i+1 . Inverting this matrix in the usual way one 
obtains a matrix relating the x values of the distribution function to the coefficients making 
it possible to compute them just from the knowledge of the different x values and the value of 
the distribution function at those x values. This calculation is implemented in NEWARRAY 
where the initial distribution is handed to the subroutine and the coefficient array is then 
returned. The coefficient array in which the values of the coefficients for the integration 
are stored, has 3 times the size of the user-specified number of points in x since we have 3 
coefficients for each bin in x. We treat the last integration between the points x and X\ 
again by approximating the distribution in this last bin through a second order polynomial. 
However, for this last bin, the coefficients are computed using the last three values in x and 
of the distribution at those points, since the point x_i which would be required according 
to the above prescription for calculating the coefficients, does not exist. 

After having regrouped the terms appearing in the convolution integral in such a way 
that all the necessary cancellations of large terms occur within the analytic expression for 
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the integral and not between different parts of the convolution integral, the integration of the 
different terms is performed in the new subroutine NINTEGR with the aid of the coefficient 
array from NEWARRAY. As mentioned above the convolution integral from x\ to 1 is split 
up into several intervals in which the integration is carried out analytically. To give an 
example of this procedure we consider the convolution integral of P qg (xi/yi, A/yi) with the 
parton distribution gs(yi)' 

f 1 dyi f l dy 1 x\{ Xl - A) r 1 d Vl Xl ( yi - Xl ) 2 

/ p qg9s = / - ( x^9s{yi) + / - ( TW9s(yi) (8) 

J xi yi Jx! yi yi(yi-A) 2 J X1 y x yi ( yi -A) 2 

suppressing presently irrelevant factors in front of the integral. The two parts in Eq. 8 
are calculated in different parts of NINTEGR and then put together in either NSRHSM, 
NSRHSP or SNRHS. In NINTEGR the integrals are split up according to Eq. 6 and then 
analytically evaluated in the different x-bins [24]. If the dependence of the integrand on A 
is only of a multiplicative nature it is enough to compute the integral for each bin once. To 
get the value of convolution integral for a term with such a A [25] dependence, it is enough 
to store the result of the integration in the bin from Xn-i to xn in the output array for this 
term at the position N — 1 [26], add to this result the value of the integral in the bin from 
xat_ 2 to xm-i and store it at the position N — 2 and so forth. In this manner one only has 
to calculate N — 1 integrals, however if the integrand has a more complicated dependence 
on A like Xi — A one needs to compute N(N — l)/2 integrals. For example in order to find 
the integration value for the xn-i bin with x\ = xjv-i one needs only one integral but at 
xn-2 we have to redo our integral for the xn-i bin since x\ = xn-2 plus we need to add the 
contribution from the xn_ 2 bin to get the correct answer for the output array at position 
iV — 2 and so forth. This need for additional evaluations of integrals slows the program down 
but in the end it turns out to be only about a factor of 4 — 5 slower ( as tested on a SPARC 
10 ) than the original CTEQ-code which is speed optimized. The integral with the regular + 
- prescription is evaluated using the routine HINTEG from the original CTEQ-code whereas 
the generalized + - prescription is evaluated according to the methods described above due 
to its nontrivial dependence on x\ and A. 



The case x\ = A = x and A << xi, are implemented in NINTEGR in the same way 
as above but separately from each other and from the more general case. For x\ = A = x 
the form of the integrands simplify in such a way that one can use the integration routines 
INTEGR and HINTEG from the original CTEQ-code. In the case of A << x\ the analytic 
expressions obtained for the above general case are expanded to first order in A and then 
the same methods as above for evaluating the integrals are applied. The last case also allows 
us to go to the diagonal case by setting A = without using the integration routines from 
the original CTEQ-code giving us a valuable tool to compare our code to the original one. 

B. Modifications in NSRHSM, NSRHSP and SNRHS 

The modification in the already existing routines NSRHSM, NSRHSP and SNRHS of 
the original CTEQ package are rather trivial. The most notable difference is that the 
subroutine NEWARRAY is called every time either of the three subroutines is called since 
the distribution function handed down on an array changes with every call of NSRHSM, 
NSRHSP and SNRHS. In NSRHSM and NSRHSP, NEWARRAY is only called once since 
one is only dealing with the non-singlet part containing no gluons, whereas in SNRHS the 
subroutine for the singlet case, one needs a coefficient array for both the quark and the 
gluon. Besides this change, the calls for INTEGR are replaced by NINTEGR according to 
how the convolution integral has been regrouped as explained in Sec. Ill A. The different 
regrouped expressions are then added, after integration for different x-values, to obtain the 
final answer in an output array which is handed back to the subroutine EVOLVE. The 
method is the same as in the original CTEQ-code but the terms themselves have changed 
of course. 

IV. CODE ANALYSIS 

As a first step we tested the stability and speed of convergence of the code and found that 
by doubling the number of points in the x-grid, which is only relevant for the convolution 
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integral, from 50 to 300 the result of our calculation changed by less than 0.5%, hence we can 
assume that our code converges rather rapidly. We also found the code to be stable down 
to an x 2 = 10 -10 beyond which we did not test. Furthermore we can reproduce the result of 
the original CTEQ-code, i.e. the diagonal case in LO within 0.5% giving us confidence that 
our code works well since the analytic expressions for the diagonal case are the expansions 
of the general case of non vanishing asymmetry up to, but not including, 0(A 2 ). 

In the following figures (Fig. 2-7) we compare, for illustrative purposes, the diagonal and 
nondiagonal case by plotting the ratio 



for various values of x±, Q 2 and A = x B j [27] , i.e. varying x 2 , using the CTEQ4M and 
CTEQ4LQ [28] parameterizations [30]. We assume the same initial conditions for the diag- 
onal and nondiagonal case (see Ref. [5] for a detailed physical motivation of this ansatz). 

The reader might wonder why only CTEQ4M and CTEQ4LQ and not GRV or MRS were 
used. The answer is not a prejudice of the authors against GRV or MRS but rather the fact 
that a comparison of CTEQ4M and CTEQ4LQ shows the same characteristic as comparing, 
for example, CTEQ4M and GRV at LO. The observation is the following: CTEQ4LQ is 
given at a different, rather low, Q, as compared to CTEQ4M and hence one has significant 
corrections from NLO terms in the evolution at these scales. This leads to a large difference 
between CTEQ4LQ and CTEQ4M (see Fig. 8), if one evolves the CTEQ4LQ set from its very 
low Q scale to the scale at which the CTEQ4M distribution is given, making a sensitivity 
study of nondiagonal parton distributions for different initial distributions impossible at 
LO. Of course, the inclusion of the NLO terms corrects this difference in the diagonal case 
but since there is no NLO calculation of the nondiagonal case available yet, a study of the 
sensitivity of nondiagonal evolution to different initial distributions has to wait. 

The figures themselves suggest the following. The lower the starting scale, the stronger 
the effect of the difference of the nondiagonal evolution as compared to the diagonal one and 



R q {xi, x 2 , Q 2 ) 



R g (x 1 ,x 2 ,Q 2 ) 



g(x 1} x 2 ,Q 2 ) 
x 1 G(x 1 , Q 2 ) 
q{xi,x 2 ,Q 2 ) 
x 1 Q(x 1 ,Q 2 ) ' 



(9) 
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also that most of the difference between nondiagonal and diagonal evolution stems from the 
first few steps in the evolution at lower scales. Secondly, under the assumption that the NLO 
evolution in the nondiagonal case will yield the same results for the parton distributions at 
some scale Q, irrespective of the starting scale Q , in analogy to the diagonal case. One 
can say that the NLO corrections to the nondiagonal evolution will be in the same direction 
and same order of magnitude as the diagonal NLO evolution. If, in the nondiagonal case, 
the NLO corrections were in the opposite direction, which would lead to a marked deviation 
from the LO results, compared to the diagonal case, the overall sign of the NLO nondiagonal 
kernels would have to change for some A ^ since in the limit A — > we have to recover 
the diagonal case. This occurance is not likely for the following reason: First, the Feynaman 
diagrams involved in the calculation of the NLO nondiagonal kernels are the same as in the 
diagonal case, except for the different kinematics, therefore, we have a very good idea about 
the type of terms appearing in the kernels, namely polynomials, logs and terms in need of 
regularization such as ln(z) ln ^~^ ■ Moreover, the kernels, as stated before, have to reduce 
to the diagonal case in the limit of vanishing A which fixes the sign of most terms in the 
kernel, thus the only type of terms which are allowed and could change the overall sign of 
the kernel are of the form 

-/(xi/i/i,A/ yi ) (10) 

which will be numerically small unless y\ ~ A in the convolution integral of the evolution 
equations. Moreover, we know that in this limit the contribution of the regularized terms 
in the kernel give the largest contributions in the convolution integral and therefore sign 
changing contributions in the nondiagonal case would have to originate from regularized 
terms. This in turn disallows a term like Eq. 10 due to the fact that regularized terms are 
not allowed to vanish in the diagonal limit, since the regularized terms arise from the same 
Feynman diagrams in the both diagonal and nondiagonal case. Therefore, the overall sign 
of the contribution of the NLO nondiagonal kernels will be the same as in the diagonal case. 
A word should be said about how the results of Ref. [10] compare to ours. For the case 
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of the same A = 1CT 3 similar starting scales and almost identical values of Q we find good 
agreement with their numbers for R g at X\ ~ A [29] and are slightly higher at larger x\. 
The observed differences are due to the fact that the quark distributions are included in 
our evolution as compared to [10] and their initial distributions is slightly different. We 
also find very similar ratios to [10] if one changes the starting scale to a lower one. The 
slight difference of a few percent in the ratios between us and [10] can again be attributed 
to the fact that they used the GRV distribution as compared to our use of the CTEQ4 
distributions, hence a slight difference in the starting scales and their lack of incorporating 
quarks into the evolution. 

V. CONCLUSIONS 

We modified the original CTEQ-code in such a way that we can now compute the evo- 
lution of nondiagonal parton distributions to LO. We gave a detailed account of the mod- 
ifications and the methods employed in the new or modified subroutines. As the reader 
can see, the modifications and methods themselves are not something magical but rather 
a straightforward application of well known numerical methods. We further demonstrated 
the rapid convergence and stability of our code. In the limit of vanishing asymmetry we 
reproduce the diagonal case in LO as obtained from the original CTEQ-code within 1%. We 
also have good agreement with the results in Ref. [10]. In the future, after the NLO kernels 
for the nondiagonal case have been calculated, we will extend the code to the NLO level to 
be on par with the diagonal case. 
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FIG. 2. R g is plotted versus x\ for fixed A using the CTEQ4M parameterization with 
=1.6 GeV and A = 202 MeV. 
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FIG. 3. R q is plotted versus x\ for fixed A using the CTEQ4M parameterization with 
=1.6 GeV and A =202 MeV. 
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FIG. 4. R g and R q are plotted versus x\ for A = using the CTEQ4M parameterization with 
=1.6 GeV and A = 202 MeV. 
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FIG. 5. R g is plotted versus x\ for fixed A using the CTEQ4LQ parameterization with 
=0.7 GeV and A =174 MeV. 



17 



4.5 



3.5 




q(x,, A,Q)/x,Q(x,,Q) at A= 1 0" 
input — CTEQ4LQ 
Q = 2.56 GeV 
Q = 5.35 GeV 
Q = 1 0.0 GeV 



10 



4.5 



q(x,, A,Q)/x,Q(x,,Q) at A= 1 0" 
input — CTEQ4LQ 
Q = 2.56 GeV 
Q = 5.55 GeV 
Q = 1 0.0 GeV 




FIG. 6. R q is plotted versus x\ for fixed A using the CTEQ4LQ parameterization with 
=0.7 GeV and A =174 MeV. 
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FIG. 7. R g and R q are plotted versus x\ for A = using the CTEQ4LQ parameterization with 
=0.7 GeV and A =174 MeV. 
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FIG. 8. The ratios for CTEQ4M to CTEQ4LQ for gluons and quarks in the diagonal case is 
plotted to demonstrate the difference between the LO evolution for these parameterizations. 
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